Set Up

Install/load necessary R packages

Set date range

min_year = 2017
max_year = 2021

Peatland Daily Water Table (Bogwell Data)

More information about peatland and upland water elevation data collected at the Marcell Experimental Forest can be found on the Environmental Data Initiative Repository Portal.

Manual Stripcharts

  • Read in the most up-to-date data from the following Box file path: External-MEF_DATA/Hydro/WaterTableElevation/Bogwells/L1

    • Less up-to-date, already published data can be found at the following Box file path: External-MEF_DATA/EDI/peatland_water_table_daily/edi.562.2/data_objects/

  • Clean the data into a tidy format

    • Subset data to 2017 - present
#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/Hydro/WaterTableElevation/Bogwells/L1"

#read in the raw data 
bogwell <- read_csv(here(filepath, "L1_Peatland_well_daily_history_2021forMia.csv"))

#clean the data
##note: we do not want any E values for the flag column 
bogwell_clean <- bogwell %>% 
  clean_names() %>% 
  mutate(year = format(as.POSIXct(date, format = '%Y-%m-%d', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year)) %>% 
  subset(year >= min_year & year <= max_year)

#save clean data CSV to use in future RMD files
write.csv(x = bogwell_clean, 
          file = file.path(here("intermediate_data", "bogwell_clean.csv")),
          row.names = FALSE)

Plot all sites together

#plot
p_bogwell <- ggplot(data = bogwell_clean) + 
  geom_line(aes(x = date, y = wte, color = peatland)) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m above sea level)",
       title = paste0("Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")"),
       color = "Site"
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7))

#static plot 
p_bogwell

S1 Manual Stripchart

bogwell_clean_s1 <- bogwell_clean %>% 
    subset(peatland == "S1")

#plot
p_bogwell_s1 <- ggplot(data = bogwell_clean_s1) + 
  geom_line(aes(x = date, 
                y = wte, 
                group = 1,
                text = paste0("Date: ", date, "\n",
                              "Water Table Elevation (m): ", wte))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m above sea level)",
       title = paste0("S1 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")") 
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7))

#interactive plot 
ggplotly(p_bogwell_s1,
         tooltip = "text")

S2 Manual Stripchart

bogwell_clean_s2 <- bogwell_clean %>% 
    subset(peatland == "S2")

#plot
p_bogwell_s2 <- ggplot(data = bogwell_clean_s2) + 
  geom_line(aes(x = date, 
                y = wte, 
                group = 1,
                text = paste0("Date: ", date, "\n",
                              "Water Table Elevation (m): ", wte))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m above sea level)",
       title = paste0("S2 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")") 
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7))

#interactive plot 
ggplotly(p_bogwell_s2,
         tooltip = "text")
bogwell_clean_s2 %>% 
  #remove rows that have NA wte values 
  filter(!is.na(wte)) %>% 
  summarise(mean_wte_m = round(mean(wte), digits = 3),
            max_wte_m = max(wte),
            min_wte_m = min(wte), 
            sd_wte_m = round(sd(wte), digits = 3)
            ) %>% 
  #t() %>% 
  kable(caption = "S2 Stripchart Bogwell Statistics", 
        col.names = c("Mean WTE (m)", 
                      "Max WTE (m)",
                      "Min WTE (m)",
                      "SD")) %>% 
  kable_material(c("striped", "hover"))
S2 Stripchart Bogwell Statistics
Mean WTE (m) Max WTE (m) Min WTE (m) SD
421.88 422.15 421.48 0.095

S3 Manual Stripchart

bogwell_clean_s3 <- bogwell_clean %>% 
    subset(peatland == "S3")

#plot
p_bogwell_s3 <- ggplot(data = bogwell_clean_s3) + 
  geom_line(aes(x = date, 
                y = wte, 
                group = 1,
                text = paste0("Date: ", date, "\n",
                              "Water Table Elevation (m): ", wte))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m above sea level)",
       title = paste0("S3 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")") 
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7))

#interactive plot 
ggplotly(p_bogwell_s3,
         tooltip = "text")

S4 Manual Stripchart

bogwell_clean_s4 <- bogwell_clean %>% 
    subset(peatland == "S4")

#plot
p_bogwell_s4 <- ggplot(data = bogwell_clean_s4) + 
  geom_line(aes(x = date, 
                y = wte, 
                group = 1,
                text = paste0("Date: ", date, "\n",
                              "Water Table Elevation (m): ", wte))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m above sea level)",
       title = paste0("S4 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")") 
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7))

#interactive plot 
ggplotly(p_bogwell_s4,
         tooltip = "text")

S5 Manual Stripchart

bogwell_clean_s5 <- bogwell_clean %>% 
    subset(peatland == "S5")

#plot
p_bogwell_s5 <- ggplot(data = bogwell_clean_s5) + 
  geom_line(aes(x = date, 
                y = wte, 
                group = 1,
                text = paste0("Date: ", date, "\n",
                              "Water Table Elevation (m): ", wte))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m above sea level)",
       title = paste0("S5 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")") 
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7))

#interactive plot 
ggplotly(p_bogwell_s5,
         tooltip = "text")

S6 Manual Stripchart

bogwell_clean_s6 <- bogwell_clean %>% 
    subset(peatland == "S6")

#plot
p_bogwell_s6 <- ggplot(data = bogwell_clean_s6) + 
  geom_line(aes(x = date, 
                y = wte, 
                group = 1,
                text = paste0("Date: ", date, "\n",
                              "Water Table Elevation (m): ", wte))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m above sea level)",
       title = paste0("S6 Daily Peatland Water Table Elevation (", min_year, "-", max_year, ")") 
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7))

#interactive plot 
ggplotly(p_bogwell_s6,
         tooltip = "text")

Shaft Encoders

Function for loading Campbell DataLogger .dat files:

readCDL = function(file){
  
  # read data file starting on 5th line
  dat <- read.csv(file, sep=",",header=FALSE,skip=4,na.strings="NAN",stringsAsFactors=F)
  
  # Read in just the header line (l2)
  # unlist the line, and remove quotes 
  h <- readLines(file, n=2)[2]
  n <- as.factor(unlist(strsplit(h, ",")) )
  n2 <- gsub('"', "", n)
  
  # assign column names to dataframe
  colnames(dat) = n2
  
  return(dat)
}

Constants to convert between feet and meters:

MetersPerFoot = 0.3048
FeetPerMeter = 3.28048

S1 Shaft Encoder

  • Read in the 2022 data from the following Box file path: External-MEF_DATA/DataDump/RealTimeData

  • Clean the 2022 data into a tidy format

    • Subset data to 2020 - present
#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/RealTimeData"

#read in the raw data 
##note: this data is only 2022 
s1_se <- readCDL(file = here(filepath, "S1-EM3_Table1.dat"))

#clean the data 
s1_se_clean <- s1_se %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         wt_elev) %>% 
  #rename columns 
  rename(datetime = timestamp,
         wt_elev_ft = wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
    wt_elev_ft = as.numeric(wt_elev_ft),
    wt_elev_m = wt_elev_ft * MetersPerFoot)

2019 data from the annual_appended_logger_files directory in Box

#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/annual_appended_logger_files"

#set the year
year = "2019"

#set the file name 
file = "S1-EM3_Table1.dat"

#read in 2019 data 
s1_2019 <- readCDL(file = here(filepath, year, file))

#clean 2019 data 
s1_2019_clean <- s1_2019 %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         wt_elev) %>% 
  #rename columns 
  rename(datetime = timestamp,
         wt_elev_ft = wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
    wt_elev_ft = as.numeric(wt_elev_ft),
    wt_elev_m = wt_elev_ft * MetersPerFoot)

2020 data from the annual_appended_logger_files directory in Box

#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/annual_appended_logger_files"

#set the year
year = "2020"

#set the file name 
file = "S1-EM3_Table1.dat"

#read in 2020 data
#note: this data is 2020 - 2022 
s1_2020 <- readCDL(file = here(filepath, year, file))

#clean 2020 data 
s1_2020_clean <- s1_2020 %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         wt_elev) %>% 
  #rename columns 
  rename(datetime = timestamp,
         wt_elev_ft = wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
    wt_elev_ft = as.numeric(wt_elev_ft),
    wt_elev_m = wt_elev_ft * MetersPerFoot)

Combine all S1 data from both the RealTimeData folder and annual_appnded_logger_files folder in Box

s1_all <- rbind(s1_se_clean, s1_2019_clean, s1_2020_clean) %>% 
  #remove NA values
  filter(!is.na(wt_elev_m)) %>% 
  #extract year 
  mutate(year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                             tz = "GMT"),
                  format = '%Y'),
    year = as.numeric(year)) %>% 
  #remove duplicate entries 
  unique()

#set min and max year 
min_year = min(s1_all$year)
max_year = max(s1_all$year)

#plot
p_s1 <- ggplot(data = s1_all) + 
  geom_line(aes(x = datetime,
                y = wt_elev_m,
                group = 1,
                text = paste0("Date: ", datetime, "\n",
                              "Water Table Elevation (m): ", round(wt_elev_m, digits = 3)))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m above sea level)",
       title = paste0("S1 Shaft Encoder Water Table Elevation (", min_year, "-", max_year, ")") 
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7))

#interactive plot 
ggplotly(p_s1,
         tooltip = "text")

S2 Shaft Encoder

  • Read in the 2022 data from the following Box file path: External-MEF_DATA/DataDump/RealTimeData
    • We are interested in the MaxWTElev column
  • Clean the 2022 data into a tidy format
  • Subset data to 2020 - present
#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/RealTimeData"

#read in the raw data 
s2_se <- readCDL(file = here(filepath, "S2_bogwell_S2BW_Daily.dat"))

#clean the data
s2_se_clean <- s2_se %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         max_wt_elev) %>% 
  #rename columns 
  rename(datetime = timestamp,
         max_wt_elev_ft = max_wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
    max_wt_elev_ft = as.numeric(max_wt_elev_ft),
    max_wt_elev_m = max_wt_elev_ft * MetersPerFoot, 
    year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year))

#set min and max year 
min_year = min(s2_se_clean$year)
max_year = max(s2_se_clean$year)
  • Plot
#plot 
p_s2 <- ggplot(data = s2_se_clean) + 
  geom_line(aes(x = datetime, 
                y = max_wt_elev_m,
                group = 1,
                text = paste0("Date: ", datetime, "\n",
                              "Water Table Elevation (m): ", max_wt_elev_m))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m)", 
       title = paste0("S2 Bogwell Water Table Elevation (", min_year, "-", max_year, ")")
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
          plot.subtitle = element_text(hjust = 0.5))

#interactive plot
ggplotly(p_s2,
         tooltip = "text")

  • Calculate S2 statistics of the maximum daily WTE including:

    • mean

    • max

    • min

    • SD

s2_se_clean %>% 
  #rename columns 
  rename(wte = max_wt_elev_m) %>% 
  #remove rows that have NA wte values 
  filter(!is.na(wte)) %>% 
  summarise(mean_wte_m = round(mean(wte), digits = 3),
            max_wte_m = max(wte),
            min_wte_m = min(wte), 
            sd_wte_m = round(sd(wte), digits = 3)
            ) %>% 
  #t() %>% 
  kable(caption = "S2 Shaft Encoder Bogwell Daily Max WTE Statistics", 
        col.names = c("Mean WTE (m)", 
                      "Max WTE (m)",
                      "Min WTE (m)",
                      "SD")) %>% 
  kable_material(c("striped", "hover"))
S2 Shaft Encoder Bogwell Daily Max WTE Statistics
Mean WTE (m) Max WTE (m) Min WTE (m) SD
421.844 422.0916 421.4854 0.114

Next, calculate the difference from the mean and SD of the difference from the mean

#create a dataframe to calculate the difference from the mean WTE (in meters)
diff <- s2_se_clean %>% 
  #rename columns 
  rename(wte = max_wt_elev_m) %>% 
  #remove rows that have NA wte values 
  filter(!is.na(wte)) %>% 
  mutate(diff_from_mean = wte - mean(wte))

#calculate the SD of the differences from the mean WTE (in meters)
sd_diff_from_mean <- sd(diff$diff_from_mean)

#set min and max year 
min_year = min(diff$year)
max_year = max(diff$year)

#plot 
ggplot(data = diff) + 
  geom_histogram(aes(x = diff_from_mean),
                 binwidth = 0.025) + 
  theme_classic() + 
  labs(x = "Differences from the Mean WTE (m)", 
       y = "Frequency", 
       title = paste0("S2 Daily Max WTE Differences from the Mean (", min_year, "-", max_year, ")")
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
          plot.subtitle = element_text(hjust = 0.5))

S3 Shaft Encoder

S3 was omitted due to unreliable data

S4 Shaft Encoder

S4_bogwell_S4BW_Daily.dat(MaxWTElev; annual_appended)

2020 data

#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/annual_appended_logger_files"

year = "2020"

file = "S4_bogwell_S4BW_Daily.dat"

#read in 2020 data 
s4_2020 <- readCDL(file = here(filepath, year, file))

#clean the data 
s4_2020_clean <- s4_2020 %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         max_wt_elev,
         max_wt_elev_m) %>% 
  #rename columns 
  rename(datetime = timestamp,
         max_wt_elev_ft = max_wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
    max_wt_elev_ft = as.numeric(max_wt_elev_ft),
    max_wt_elev_m = as.numeric(max_wt_elev_m), 
    date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year))

#set min and max year 
min_year = min(s4_2020_clean$year)
max_year = max(s4_2020_clean$year)
  • Plot
#plot 
p_s4 <- ggplot(data = s4_2020_clean) + 
  geom_line(aes(x = datetime, 
                y = max_wt_elev_m,
                group = 1,
                text = paste0("Date: ", date, "\n",
                              "Water Table Elevation (m): ", max_wt_elev_m))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m)", 
       title = paste0("S4 Bogwell Water Table Elevation (", min_year, "-", max_year, ")")
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
          plot.subtitle = element_text(hjust = 0.5))

#interactive plot
ggplotly(p_s4,
         tooltip = "text")

S5 Shaft Encoder

2020-2022 data

#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/annual_appended_logger_files"

year = "2020"

file = "S5_bogwell_S5BW_Daily.dat"

#read in 2020-2022 data 
s5_2020 <- readCDL(file = here(filepath, year, file))

#clean the data 
s5_2020_clean <- s5_2020 %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         max_wt_elev,
         max_wt_elev_m) %>% 
  #rename columns 
  rename(datetime = timestamp,
         max_wt_elev_ft = max_wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
    max_wt_elev_ft = as.numeric(max_wt_elev_ft),
    max_wt_elev_m = as.numeric(max_wt_elev_m), 
    date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year))

RealTimeData

#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/RealTimeData"

#read in the raw data 
##note: this data is 2021 - 2022 
s5_se <- readCDL(file = here(filepath, "S5_bogwell_S5BW_Daily.dat"))

#clean the data 
s5_se_clean <- s5_se %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         max_wt_elev,
         max_wt_elev_m) %>% 
  #rename columns 
  rename(datetime = timestamp,
         max_wt_elev_ft = max_wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
    max_wt_elev_ft = as.numeric(max_wt_elev_ft),
    max_wt_elev_m = as.numeric(max_wt_elev_m), 
    date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year))

Combine all S5 data and plot

s5_all <- rbind(s5_2020_clean, s5_se_clean) %>% 
  #remove NA values
  filter(!is.na(max_wt_elev_m)) %>% 
  #remove duplicate observations 
  distinct()

#set min and max year 
min_year = min(s5_se_clean$year)
max_year = max(s5_se_clean$year)

#plot 
p_s5 <- ggplot(data = s5_all) + 
  geom_line(aes(x = datetime, 
                y = max_wt_elev_m,
                group = 1,
                text = paste0("Date: ", date, "\n",
                              "Water Table Elevation (m): ", max_wt_elev_m))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m)", 
       title = paste0("S5 Bogwell Water Table Elevation (", min_year, "-", max_year, ")")
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
          plot.subtitle = element_text(hjust = 0.5))

#interactive plot
ggplotly(p_s5,
         tooltip = "text")

S6 Shaft Encoder

RealTimeData

#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/RealTimeData"

#read in the raw data 
##note: this data is 2021 - 2022 
s6_se <- readCDL(file = here(filepath, "S6_bogwell_S6BW_Daily.dat"))

#clean the data 
s6_se_clean <- s6_se %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         max_wt_elev,
         max_wt_elev_m) %>% 
  #rename columns 
  rename(datetime = timestamp,
         max_wt_elev_ft = max_wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
    max_wt_elev_ft = as.numeric(max_wt_elev_ft),
    max_wt_elev_m = as.numeric(max_wt_elev_m), 
    date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year))

2019 data

#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/annual_appended_logger_files"

year = "2019"

file = "S6_bogwell_S6BW_Daily.dat"

#read in 2019 data 
s6_2019 <- readCDL(file = here(filepath, year, file))

#clean 2019 data 
s6_2019_clean <- s6_2019 %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         max_wt_elev,
         max_wt_elev_m) %>% 
  #rename columns 
  rename(datetime = timestamp,
         max_wt_elev_ft = max_wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
         date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year),
    max_wt_elev_ft = as.numeric(max_wt_elev_ft),
    max_wt_elev_m = as.numeric(max_wt_elev_m))

2020 data

#set file path to read in data from Box 
filepath = "/Users/miaforsline/Library/CloudStorage/Box-Box/External-MEF_DATA/DataDump/annual_appended_logger_files"

year = "2020"

file = "S6_bogwell_S6BW_Daily.dat"

#read in 2020-2022 data 
s6_2020 <- readCDL(file = here(filepath, year, file))

#clean 2020-2022 data 
s6_2020_clean <- s6_2020 %>% 
  clean_names() %>% 
  #select relevant columns
  select(timestamp,
         max_wt_elev,
         max_wt_elev_m) %>% 
  #rename columns 
  rename(datetime = timestamp,
         max_wt_elev_ft = max_wt_elev) %>% 
  #reformat columns into appropriate classes 
  mutate(datetime = as.POSIXct(datetime, format = '%Y-%m-%d %H:%M:%S', 
                                  tz = "GMT"),
         date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year),
    max_wt_elev_ft = as.numeric(max_wt_elev_ft),
    max_wt_elev_m = as.numeric(max_wt_elev_m))

Combine S6 data

s6_all <- rbind(s6_se_clean, s6_2019_clean, s6_2020_clean) %>% 
  #remove one row where WTE = 0 
  subset(max_wt_elev_m > 0) %>% 
  #remove NA values
  filter(!is.na(max_wt_elev_m))

#set min and max year 
min_year = min(s6_all$year)
max_year = max(s6_all$year)

#plot
p_s6 <- ggplot(data = s6_all) + 
  geom_line(aes(x = datetime,
                y = max_wt_elev_m,
                group = 1,
                text = paste0("Date: ", date, "\n",
                              "Max Water Table Elevation (m): ", round(max_wt_elev_m, digits = 3)))) + 
  theme_classic() + 
  labs(x = "Time", 
       y = "Water Table Elevation (m above sea level)",
       title = paste0("S6 Shaft Encoder Water Table Elevation (", min_year, "-", max_year, ")") 
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7))

#interactive plot 
ggplotly(p_s6,
         tooltip = "text")

1:1 Plots

S1 1:1 Plot

Format data

#format stripchart data
bog1 <- bogwell_clean_s1 %>% 
  #remove unnecessary columns 
  select(-peatland, -flag) %>% 
  #create a column to note what type of data this is 
  mutate(type = "stripchart")

#format shaft encoder data 
se1 <- s1_all %>% 
  #extract date and year 
  mutate(date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year),
         #note data type
         type = "shaft_encoder") %>% 
  #rename columns to join the data  
  rename(wte = wt_elev_m) %>% 
  group_by(date) %>% 
  #calculate max WTE per day 
  mutate(max_wte = max(wte)) %>% 
  #remove unnecessary columns 
  select(-datetime, -wt_elev_ft, -wte) %>% 
  distinct() %>% 
  #rename columns to join the data
  rename(wte = max_wte)

#join the data 
s1_join <- full_join(x = bog1, 
                     y = se1,
                     by = c("wte", "date", "year", "type"))

#subset the join to the time period where the two data sources overlap 
s1_join_sub <- s1_join %>% 
  pivot_wider(names_from = "type",
              values_from = "wte") %>% 
  #remove rows with NA values 
  filter(!is.na(shaft_encoder), 
         !is.na(stripchart))

#set min and max year 
min_year = min(s1_join_sub$year)
max_year = max(s1_join_sub$year)

Plot

p1 <- ggplot(data = s1_join_sub) + 
  geom_point(aes(x = stripchart, 
                 y = shaft_encoder,
                 alpha = 0.1,
                 group = 1,
                 text = paste0("Date: ", date, "\n",
                               "Stripchart: ", stripchart, "\n",
                              "Shaft Encoder: ", round(shaft_encoder, digits = 3)))) + 
  #plot straight line to compare point spread 
  geom_abline(intercept = 0,
              slope = 1) + 
  theme_classic() + 
  labs(x = "Manual Stripchart", 
       y = "Shaft Encoder",
       title = paste0("S1 1:1 Daily Max WTE Comparison (", min_year, "-", max_year, ")")
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7),
        legend.position = "none")

#interactive plot 
ggplotly(p = p1, 
         tooltip = "text")

S2 1:1 Plot

Format data

#format stripchart data
bog <- bogwell_clean_s2 %>% 
  #remove unnecessary columns 
  select(-peatland, -flag) %>% 
  #create a column to note what type of data this is 
  mutate(type = "stripchart")

#format shaft encoder data 
se <- s2_se_clean %>% 
  #extract date and year 
  mutate(date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year),
         #note data type
         type = "shaft_encoder") %>% 
  #remove unnecessary columns 
  select(-datetime, -max_wt_elev_ft) %>% 
  #rename columns to join the data  
  rename(wte = max_wt_elev_m)

#join the data 
s2_join <- full_join(x = bog, 
                     y = se,
                     by = c("wte", "date", "year", "type"))

#subset the join to the time period where the two data sources overlap 
s2_join_sub <- s2_join %>% 
  pivot_wider(names_from = "type",
              values_from = "wte") %>% 
  #remove rows with NA values 
  filter(!is.na(shaft_encoder), 
         !is.na(stripchart))

#set min and max year 
min_year = min(s2_join_sub$year)
max_year = max(s2_join_sub$year)

Plot

p2 <- ggplot(data = s2_join_sub) + 
  geom_point(aes(x = stripchart, 
                 y = shaft_encoder,
                 alpha = 0.1,
                 group = 1,
                 text = paste0("Date: ", date, "\n", 
                               "Stripchart: ", stripchart, "\n",
                              "Shaft Encoder: ", round(shaft_encoder, digits = 3)))) + 
  #plot straight line to compare point spread 
  geom_abline(intercept = 0,
              slope = 1) + 
  theme_classic() + 
  labs(x = "Manual Stripchart", 
       y = "Shaft Encoder",
       title = paste0("S2 1:1 WTE Comparison (", min_year, "-", max_year, ")")
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7),
        legend.position = "none")

#interactive plot 
ggplotly(p = p2, 
         tooltip = "text")

S4 1:1 Plot

Format data

#format stripchart data
bog4 <- bogwell_clean_s4 %>% 
  #remove unnecessary columns 
  select(-peatland, -flag) %>% 
  #create a column to note what type of data this is 
  mutate(type = "stripchart")

#format shaft encoder data 
se4 <- s4_2020_clean %>% 
  #extract date and year 
  mutate(date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year),
         #note data type
         type = "shaft_encoder") %>% 
  #remove unnecessary columns 
  select(-datetime, -max_wt_elev_ft) %>% 
  #rename columns to join the data  
  rename(wte = max_wt_elev_m)

#join the data 
s4_join <- full_join(x = bog4, 
                     y = se4,
                     by = c("wte", "date", "year", "type"))

#subset the join to the time period where the two data sources overlap 
s4_join_sub <- s4_join %>% 
  pivot_wider(names_from = "type",
              values_from = "wte") %>% 
  #remove rows with NA values 
  filter(!is.na(shaft_encoder), 
         !is.na(stripchart))

#set min and max year 
min_year = min(s4_join_sub$year)
max_year = max(s4_join_sub$year)

Plot

p4 <- ggplot(data = s4_join_sub) + 
  geom_point(aes(x = stripchart, 
                 y = shaft_encoder,
                 alpha = 0.1,
                 group = 1,
                 text = paste0("Date: ", date, "\n", 
                               "Stripchart: ", stripchart, "\n",
                              "Shaft Encoder: ", round(shaft_encoder, digits = 3)))) + 
  #plot straight line to compare point spread 
  geom_abline(intercept = 0,
              slope = 1) + 
  theme_classic() + 
  labs(x = "Manual Stripchart", 
       y = "Shaft Encoder",
       title = paste0("S4 1:1 WTE Comparison (", min_year, "-", max_year, ")")
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7),
        legend.position = "none")

#interactive plot 
ggplotly(p = p4, 
         tooltip = "text")

S5 1:1 Plot

Format data

#format stripchart data
bog5 <- bogwell_clean_s5 %>% 
  #remove unnecessary columns 
  select(-peatland, -flag) %>% 
  #create a column to note what type of data this is 
  mutate(type = "stripchart")

#format shaft encoder data 
se5 <- s5_all %>% 
  #extract date and year 
  mutate(date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year),
         #note data type
         type = "shaft_encoder") %>% 
  #remove unnecessary columns 
  select(-datetime, -max_wt_elev_ft) %>% 
  #rename columns to join the data  
  rename(wte = max_wt_elev_m)

#join the data 
s5_join <- full_join(x = bog5, 
                     y = se5,
                     by = c("wte", "date", "year", "type"))

#subset the join to the time period where the two data sources overlap 
s5_join_sub <- s5_join %>% 
  pivot_wider(names_from = "type",
              values_from = "wte") %>% 
  #remove rows with NA values 
  filter(!is.na(shaft_encoder), 
         !is.na(stripchart))

#set min and max year 
min_year = min(s5_join_sub$year)
max_year = max(s5_join_sub$year)

Plot

p5 <- ggplot(data = s5_join_sub) + 
  geom_point(aes(x = stripchart, 
                 y = shaft_encoder,
                 alpha = 0.1,
                 group = 1,
                 text = paste0("Date: ", date, "\n", 
                               "Stripchart: ", stripchart, "\n",
                              "Shaft Encoder: ", round(shaft_encoder, digits = 3)))) + 
  #plot straight line to compare point spread 
  geom_abline(intercept = 0,
              slope = 1) + 
  theme_classic() + 
  labs(x = "Manual Stripchart", 
       y = "Shaft Encoder",
       title = paste0("S5 1:1 WTE Comparison (", min_year, "-", max_year, ")")
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7),
        legend.position = "none")

#interactive plot 
ggplotly(p = p5, 
         tooltip = "text")

S6 1:1 Plot

Format data

#format stripchart data
bog6 <- bogwell_clean_s6 %>% 
  #remove unnecessary columns 
  select(-peatland, -flag) %>% 
  #create a column to note what type of data this is 
  mutate(type = "stripchart")

#format shaft encoder data 
se6 <- s6_all %>% 
  #extract date and year 
  mutate(date = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y-%m-%d'),
         date = as.POSIXct(date, tz = "GMT"),
         year = format(as.POSIXct(datetime, format = '%m/%d/%Y %H:%M:%S', 
                                  tz = "GMT"), 
                       format = '%Y'),
         year = as.numeric(year),
         #note data type
         type = "shaft_encoder") %>% 
  #remove unnecessary columns 
  select(-datetime, -max_wt_elev_ft) %>% 
  #rename columns to join the data  
  rename(wte = max_wt_elev_m)

#join the data 
s6_join <- full_join(x = bog6, 
                     y = se6,
                     by = c("wte", "date", "year", "type"))

#subset the join to the time period where the two data sources overlap 
s5_join_sub <- s5_join %>% 
  pivot_wider(names_from = "type",
              values_from = "wte") %>% 
  #remove rows with NA values 
  filter(!is.na(shaft_encoder), 
         !is.na(stripchart))

#set min and max year 
min_year = min(s5_join_sub$year)
max_year = max(s5_join_sub$year)

Plot

p5 <- ggplot(data = s5_join_sub) + 
  geom_point(aes(x = stripchart, 
                 y = shaft_encoder,
                 alpha = 0.1,
                 group = 1,
                 text = paste0("Date: ", date, "\n", 
                               "Stripchart: ", stripchart, "\n",
                              "Shaft Encoder: ", round(shaft_encoder, digits = 3)))) + 
  #plot straight line to compare point spread 
  geom_abline(intercept = 0,
              slope = 1) + 
  theme_classic() + 
  labs(x = "Manual Stripchart", 
       y = "Shaft Encoder",
       title = paste0("S5 1:1 WTE Comparison (", min_year, "-", max_year, ")")
       ) + 
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5, 
                                     size = 7),
        legend.position = "none")

#interactive plot 
ggplotly(p = p5, 
         tooltip = "text")

Note that this file is being knit as index.html into the climate_bogwell repository in order to update the GitHub pages website, where the plots can be viewed online.